今回は前々回と同じデータを生成して分析していく。
個人の所得\(Y\)と学歴\(X\)・能力\(A\)との関係についての架空データを次のように生成する(ただし\(\varepsilon\)は測定誤差などを反映した撹乱項)
\[ Y = 200 + 10A + 500X+ \varepsilon \]
# 事前準備 --------------------
# パッケージの読み込み
library(tidyverse)
# 乱数の種を固定
set.seed(0)
# データの生成 ----------------
n = 100000
# 能力は0から100まで均等に分布
ability = runif(n, min = 0, max = 100)
# IDとabilityをデータフレームに格納する
df = data_frame(ID = 1:n, ability)
# 大卒フラグ
## 条件1:能力が 80 以上の約 2 万人の中から約 1 万人がランダムに選ばれて大卒となる。
college = df %>% filter(ability >= 80) %>% sample_frac(0.5) # 大卒の人
college["is_college1"] = 1
no_college = anti_join(df, college, by = c("ID","ability")) # 大卒じゃない人
no_college["is_college1"] = 0
df = bind_rows(college, no_college) %>% arrange(ID) # 両者を結合
## 条件2:能力を部分的に反映した学力テストの点数が 180 点以上であれば大卒となる
### 1万人くらいが該当するようにする
df["score"] = 30 * log10(ability) + rnorm(n, mean = 115, sd = 10)
df["is_college2"] = 1*(df["score"] >= 180)
## 2つの条件のいずれかに該当していれば大卒とする
df = df %>% mutate(is_college = case_when(is_college1 == 1 | is_college2 == 1 ~ 1, # "|"はorの記号
TRUE ~ 0)) # それ以外は0
# 所得
df["income"] = 200 + 10*df["ability"] + 500*df["is_college"] + rnorm(n, sd = 50)
# 最初の6行
head(df)## # A tibble: 6 x 7
## ID ability is_college1 score is_college2 is_college income
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 89.7 1 172. 0 1 1572.
## 2 2 26.6 0 143. 0 0 479.
## 3 3 37.2 0 158. 0 0 492.
## 4 4 57.3 0 172. 0 0 720.
## 5 5 90.8 1 184. 1 1 1738.
## 6 6 20.2 0 170. 0 0 328.
こうして生成したデータの所得と能力・学歴・テストの点数の関係をプロットすると次のようになる。
以下では「学歴(大卒になること)が所得をどれだけ上昇させるのか」を知ることを目的として分析していくことにする。
temp <- df %>%
mutate(education = case_when(is_college2 == 1 ~ "大卒(条件2)",
is_college == 0 ~ "非大卒")) %>%
filter(!is.na(education))
ggplot(temp,
aes(x = score, y = income, color = education, fill = education))+
geom_point(alpha = 0.5)+
geom_smooth(method="lm",se=F, color="black", alpha = 0.5)+
labs(title = "大卒(条件2)と非大卒のスコアと所得")条件2で大卒になった人の場合を見ていくと,テストの点数がある閾値(180点)を超える点で所得にジャンプがみられる。 回帰不連続デザイン(regression discontinuity design: RDD)はこうした閾値を伴うDGPを利用する調査設計(デザイン)で,以下のようなロジックをもとに分析を行う。
RDDは閾値周辺のデータを使った回帰分析である局所線形回帰(local linear regression)などを使う。
# 条件2の大卒者と非大卒者のみを取り出す
df_rdd = df %>% filter(is_college2 == 1 | is_college == 0)
# 閾値周辺のサンプルの取り出し
h = 1 # バンド幅:閾値周辺のどの程度の幅のデータを使うか
c = 180 # 閾値
df_local = df_rdd %>% filter(c - h <= score & score < c + h)取り出した閾値周辺のデータの分布は次のようになっている
# plot
ggplot(df_local, aes(x = score, y = income))+
geom_point()# 局所線形回帰
llr <- lm(income ~ score + is_college, data = df_local)
stargazer::stargazer(llr, type = "text")##
## ===============================================
## Dependent variable:
## ---------------------------
## income
## -----------------------------------------------
## score 13.734
## (12.486)
##
## is_college 534.456***
## (14.439)
##
## Constant -1,581.780
## (2,241.268)
##
## -----------------------------------------------
## Observations 2,754
## R2 0.674
## Adjusted R2 0.673
## Residual Std. Error 190.169 (df = 2751)
## F Statistic 2,838.943*** (df = 2; 2751)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
比較的500に近い値を推定することができた